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Abstract 



The random field Ising model with Gaussian disorder is studied using a 
new Monte Carlo algorithm. The algorithm combines the advantanges of 
the replica exchange method and the two-replica cluster method and is much 
more efficient than the Metropolis algorithm for some disorder realizations. 
Three-dimensional sytems of size 24^ are studied. Each realization of disorder 
is simulated at a value of temperature and uniform field that is adjusted to 
the phase transition region for that disorder realization. Energy and magne- 
tization distributions show large variations from one realization of disorder 
to another. For some realizations of disorder there are three well separated 
peaks in the magnetization distribution and two well separated peaks in the 
energy distribution suggesting a first-order transition. 
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I. INTRODUCTION 



Despite twenty five years of experimental and tfieoretical effort, pliase transitions in 
systems witfi quencfied random fields are still poorly understood. The simplest theoretical 
model is the random field Ising model (RFIM). The RFIM phase transition is believed to 
be in the same universality class as the phase transitions in diluted antiferromagnets in 
a uniform field and fluids in porous media. The three dimensional RFIM is known 
to have an ordered phase at sufficiently low temperature and for weak random fields. As 
the temperature or the strength of the randomness is increased, there is a transition to a 
disordered phase. The nature of this transition is not well understood. 

In this paper we describe a new replica exchange algorithm for simulating the RFIM 
and present numerical results for systems of L x L x L spins with L up to 24. We show 
that the qualitative features of the transition differ strongly from realization to realization of 
disorder. Our results can be interpreted as suggesting that the RFIM transition is first-order 
or that a modified version of the droplet picture holds. 

The RFIM is described by the energy 

- u/kBT = (3Y. S^S, + h^S, + HY,S, (1.1) 

where the spin variables Si = ±1 reside on a lattice, /3 = l/ksT is the inverse temperature, 
the first sum is over nearest neighbor pairs on the lattice, H is an external field and A is 
the strength of the disorder. The random fields are independent random variables chosen 
from a distribution with mean zero and variance one. In this study, the random fields are 
Gaussian and the lattice is simple cubic. 

Currently, it is not known whether the phase transition for the 3D RFIM is first-order 
or continuous. Experiments on magnetic systems have been plagued by problems of poor 
equilibration and have yielded confusing results, but it appears that there is no latent heat at 
the transition and this has usually been interpreted as evidence for a continuous transition. 
Theoretical analyses |Q have also favored a continuous transition, although in many cases 
this is an assumption rather than a conclusion, and some recent work |^ suggests a fluctu- 
ation driven first-order transition. The standard picture is that of a continuous transition 
controlled by a zero temperature fixed point. The scaling theory of this transition has 
three independent critical exponents and modified hyperscaling relations. Because of the 
fixed point governing the transition has strong disorder, controlled renormalization group 
calculations have not been possible. Migdal-Kadanoff renormalization group calculations 
indicate a continuous transition PJlO| but also mistakenly predict that the g-state Potts 
transition is continuous for all q when, for 3D it is known to be first-order for g > 3. Se- 
ries analyses initially supported a first-order transition |]ll| but more recently point to a 
continuous transition, at least for weak disorder WM. Alternatively, it may be that the 



transition is continuous for weak disorder and then becomes first-order for strong disorder 
with a tricritical point separating the critical line from the first-order line. 



Recent Monte Carlo simulations ||T3|-[T5[1 have also been interpreted as showing a contin- 
uous transitions but with a jump in the magnetization. These simulations have been limited 
to system size 16^. The jump in the magnetization can be interpreted as a very small value 
of the magnetization exponent but might also signal a first-order transition. Simulations of 



2 



systems up to 64^ |T6| were interpreted as indicating a first-order transition but were clearly 



not equilibrated in the transition region. A number of numerical experiments have also been 
carried out on the RFIM at zero temperature using polynomial time ground state algorithms. 



Originally these supported a continuous transition with a small magnetization exponent |T7 



but more recent studies on larger systems |JT8|] show a jump in the magnetization and are 
thus suggestive of a first-order transition. 

The numerical results presented in this paper differ from those of previous studies in two 
significant respects. First we use an efficient algorithm that permits us to reach equilibrium 
for larger systems than those of past studies. Second, we fine-tune both the temperature and 
the external field for each realization of disorder to be as close as possible to the transition 
for that realization. 

To understand the motivation for this fine-tuning, let us suppose for the moment that the 
transition is first-order. If this is the case, then for periodic or helical boundary conditions 
we expect that there will be three phases in coexistence at the transition point. We call the 
coexisting phases +, — and 0. The ordered phases, -|- and — , have long range order and 
finite magnetization. The disordered, phase has no long range order, no magnetization 
and is characterized by spins that are predominately aligned with their local fields. The 
bond energy, defined by the first term on the RHS of Eq. ( |1.1| ), is greater in the phase 
than the ordered phases. The expected phase diagram in the T-H plane in the vicinity of 
the point of three phase coexistence is shown in Fig. ^ For T < T^. and H = there is 
phase coexistence of the + and — phases that ends at the thermal first-order transition at 
T = Tc and H = (the black dot in the figure). Since the and +(— ) phases differ in both 
energy and magnetization, the disordered phase can be maintained in coexistence with the 
-|-(— ) phase by increasing the temperature and increasing (decreasing) the field. The two 
order-disorder lines corresponding to 0+ and 0— coexistence form the arms of the "Y" in 
the figure. 

Since the disorder is independent and homogenous the free energy is self-averaging, 
hence, in the thermodynamic limit the transition occurs at a definite point (Tc, H = 0) for 
almost all realizations of the random fields. However, for finite systems, the location and 
qualitative features of the transition will depend on the realization of disorder. In particular, 
for a system of size L, with realization of random field, {h^} and disorder strength A there 
may be a point where three phases coexist at (Tc(A; {hi}), Hc{A; {hi})). A first guess, based 
on the net field due to the random fields, is that Tc(A; {hi}) and Hc{A; {hi})) should be 
displaced from the average value Tc{L) and if = by an amount of order AL'^/^. 

How accurately must H and T be fine-tuned to see three phase coexistence if it exists? 
Suppose first that we want the + and — phases to coexist and that the magnetizations of 
these phases at the transition differs by 2m. If the external field deviates from the correct 
value by 6H, the free energy difference between the + and — phase is 2m6HL'^. For both 
phases to have a significant probability the free energy difference must not greatly exceed kT. 
Thus H must be set to Hc{A; {h^}) to within an accuracy of kT/2mL'^ to have both phases 
represented. H(.{A] {h^}) will itself fluctuate from sample to sample as A/L'^/^. Similarly, 
if the entropy difference per spin between the the ordered and disordered phases is s, then 
T must be fine-tuned to within kT/sL'^ to allow to these phases to coexist. Presumably 
sample to sample fluctuations in Tc(A; {hi}) also scale as l/L'^/^. Thus if a single value of 
H = and T = Tc{L) is chosen for all realizations of the random field, one will almost never 
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see more than one phase at a time. 



II. NUMERICAL METHODS 

We use an algorithm that combines the rephca exchange method, first introduced 



by Swendsen and Wang [19| and the two- rephca cluster method of Redner, Machta and 



Chayes pO|,^. Our method is also closely related to simulated and parallel tempering p2[ 



The idea of this approach is to simultaneously simulate K replicas of the system. All repli- 
cas have the same normalized random field {hi} but each replica has different values of the 
other parameters. Replica k {k = 0, . . . , K — 1) has inverse temperature (3^, strength of 
randomness and external field H^. The replicas form a sequence so that neighboring 
replicas in the sequence are nearby in the {f3,A,H) parameter space. Neighboring replicas 
exchange magnetization with one another according to a procedure described below. One 
end of the sequence of replicas is at a value of [3, A and H that can be efficiently simulated 
using a known method while the other end of the sequence is at a value of the parameters 
that we would like to study. In our case, the replicas lie along the RFIM phase transition 
line starting from the pure Ising values /5o ~ .22615 and Aq = Hq = as shown in Fig. |^. 
The replicas are equally spaced in A. The pure Ising replica is simulated using the Wolff 
single cluster algorithm [ P5| , P^ . (We have also experimented with replicas lying along a 
line of constant (3H starting at the RFIM phase boundary, extending into the paramagnetic 
phase, and ending at a temperature high enough that the model can be efficiently simulated 
using the ordinary single-spin-fiip Metropolis algorithm. This approach, however, is found 
to be less efficient than the one described above.) 

Magnetization is exchanged between neighboring replicas using a generalization of the 
two-replica cluster method. Suppose we have two replicas at A, H) and (/?', A', H') with 
spin configurations {Si} and {S^}, respectively. A site j is considered active for this pair 
of replicas if 5*^ S'j. A bond ij between neighboring sites i and j is satisfied if Si = Sj 
and 5*^' = Sj. A cluster of active sites is formed starting from a randomly chosen active site. 
New active sites are added to the cluster by occupying satisfied bonds on the perimeter of 
the cluster with probability p{P, /?') where 

= l-e-2(^+'3'). (2.1) 

If a bond connecting a site to a cluster is occupied, the site is added to the cluster and the 
set of bonds on the perimeter is updated. In this way, the cluster grows until no further 
sites are added. The procedure is very similar to the way clusters are grown in the Wolff 
single cluster algorithm. 

Once a cluster is identified it is flipped with a probability that depends on the change 
in boundary and field energy so as to satisfy detailed balance. Flipping a cluster means 
changing the sign of all the spins in the cluster or, equivalently, exchanging the values of the 
spins in the cluster between the two replicas. The probability to fiip a cluster C with | C \ 
sites depends on the quantity S 

E = -[(A - A')hc + {H-H')\C\ +(/? - /5')(iV++ - iV__)]S'c (2.2) 

where A^++ and are, respectively, the number of -|--|- and sites that are nearest 

neighbors of the cluster, he is the net random field acting on the cluster. 
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he = 

iec 



(2.3) 



and Sc is the spin value of the cluster in the unprimed replica. If S < then the cluster is 
flipped, otherwise it is flipped with probability e~^. 

It is straightforward but tedious to show that the choice of p{P, /?') (Eq. P?T| ) together with 
the flipping probability defined by E is precisely what is needed to ensure detailed balance. 
The motivation for the choice of p{f3, j3') is most easily understand by considering the limit 
where two replicas are at the same values of /3, A and H. In that case p(/5, (3) = 1 — e~^^ and 
E = 0. Since E = 0, clusters are always flipped just as is the case for the Wolff single cluster 
algorithm. Furthermore, we have shown in Ref. pO|,^ that the active clusters percolate at 
the RFIM phase transition. If the transition is continuous clusters of all sizes are flipped. 
If the transition is flrst-order, there will be two distinct kinds of clusters; some clusters 
will be extensive and change the phase of the system while other clusters will have sizes 
less than or equal to the correlation length. In either case, the clusters identifled by the 
two-replica procedure correspond to the fluctuations that actually occur in the system at 
the phase transition and permits large changes in the spin conflguration in a single Monte 
Carlo sweep. When the two replicas do not have equal values of the parameters, then the 
clusters do not flip freely but if the replicas are close together in the parameter space, the 
acceptance fraction for flipping clusters will remain high. 

Our method and the original replica exchange method |]l9l on which it is based is similar 
to parallel tempering ||2^. In all these methods, groups of spins are exchanged between 
neighboring replicas along a sequence. In parallel tempering the whole spin conflguration 
is exchanged and the Boltzmann factor controlling the acceptance of the move depends 
on the energy difference between the replicas. In our algorithm, only some of the spins 
are exchanged and for a given distance in the parameter space between the replicas the 
acceptance fraction is larger than for parallel tempering. The consequence is that fewer, less 
closely spaced replicas are needed for the replica exchange method. 

In order to flnd the phase transition temperature and external fleld (/3c(A; {hi}), 
Hc{A; {hi})) for a given realization and strength of disorder we use a feedback mechanism. 
Starting from an initial value of (3 and H we monitor the magnetization of the system after 
each Monte Carlo sweep. If the absolute value of the magnetization is less than a lower 
cut-off Mic the system is interpreted to be in a high temperature phase and the inverse 
temperature is increased by a small amount e/3. If the absolute value of the magnetization 
is greater than an upper cut-off Muc the system is interpreted to be in one of the ordered 
phases and P is decreased by e^j. In this case, the external fleld is also adjusted by an amount 
en if the magnetization is positive and by —en if the magnetization is negative. The average 
value of P and H is computed for the period when the feedback procedure is on and is taken 
to be (/3c(A; {hi}), Hc{A; {hi})). The feedback procedure is then turned off and these values 
are used for a long equilibrium simulation. 

Most of our simulations were for 24^ systems. Except as otherwise noted, the simulations 
used K = 16 replicas equally spaced in A with the most disordered replica having A = .35. 
Initially the replicas lie on an elliptical curve in the temperature-disorder plane that starts 
at the pure Ising transition {j3 = 0.22165) and ends at the zero temperature transition 



(A//3 = 2.35) as shown in Fig. It was shown in Ref. |jT5| that this curve is a good 



estimate for the phase boundary. The feedback parameters were Mi^ = 2000, M^c = 4000, 
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e/3 = 10~^ and eu = 10~^. The feedback mechanism was run for 2 x 10^ Monte Carlo sweeps 
and then, with the temperature and external field determined by the feedback procedure, 
data were collected for 8 x 10^ Monte Carlo sweeps. Each Monte Carlo sweep is defined 
as 500 cluster steps; since the average cluster size was found to be roughly 1000, a single 
Monte Carlo sweep corresponds to attempting to flip every spin in every replica about once. 
For each cluster step a random integer k was chosen between and K. If = 0, a Wolff 
move was performed on the pure replica, if 1 < A; < fT, a two replica cluster move was 
carried out between replica k and replica k — 1. If k = K, a. two replica cluster move was 
carried out between rephcas K — 1 and K — 2. The procedure described above, totaling 
10® sweeps, was performed for 22 realizations of disorder, a number chosen in advance of 
the experiment, labeled 10 through 31 by the seed of the random number generator that 
produced {hi}. Each realization required about six days on a 450MHz Pentium III machine. 
A number of additional simulations were carried out for specific realizations of disorder (not 
all in the range 10 through 31) using parameters that differed in some way from the above. 
In particular, we did a careful study of realization 1 that is described below. 

III. RESULTS 

A. Energy and magnetization histograms 

The main lesson of our work is that each realization of disorder has its own character. 
This is best seen by examining the probability distributions of magnetization and energy 
for several realizations of disorder. Figures ^, ^ ^ and show the magnetization his- 
togram for realizations 20, 21, 25, 31 and 14, respectively. The magnetization histogram for 
realization 20, characterized by two broad maxima roughly symmetrical about the origin, is 
typical of the majority of realizations. Apart from the asymmetry in the two peaks, realiza- 
tion 20 and those like it are similar to the pure Ising model whose magnetization histogram 
is shown in Fig. |^ at the infinite system size critical temperature, /5c(0) = .22165. The 
magnetization histogram for realization 21, Fig. ^ also displays two peaks but now one peak 
is much sharper than the other and quite asymmetric about the origin. The magnetization 
histogram for realizations 25, Fig. |^ and 31, Fig. ^ are qualitatively different, displaying 
three distinct and well separated peaks. It should be noted that the feedback procedure was 
initially unsuccessful for these realizations and had to be run again with Mic = Muc = 5000 
to find all three peaks. Failure of the feedback scheme was evidenced by the presence of only 
one narrow peak in the magnetization histogram. In the case of realization 31, the feedback 
procedure with the revised parameters failed on the most disordered replica so Fig. ^ shows 
the second to last replica in the sequence at A = .3267. For both realization 25 and 31, 
the original run showed three peaked structures at weaker disorder giving a strong hint that 
the same qualitative feature would be found at stronger disorder. Finally, realization 14, 
Fig. was unique, displaying a five peaked structure. Out of a total of 22 magnetization 
histograms, 19 show two peaks, 2 three peaks and 1 five peaks. It is possible but unlikely 
that some of the two peak realizations would display more peaks with different fine-tuning. 
For example, we used other alternate values of the fine-tuning parameters to search for an 
additional positive peak for realization 21 but none was found. 
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The bond energy distribution displayed two possible qualitative behaviors. Figures 



l0| , |TT], |l^ and |l^ show the bond energy histogram for realizations 20, 21, 25, 31 and 14, 
respectively. For comparison. Fig. [1^ shows the energy distribution for the pure Ising model 
at the infinite system size critical temperature. For realization 20 and those like it, the bond 
energy distribution has a single peak and is qualitatively like that of the pure Ising model at 
criticality, except that the peak is shifted to significantly lower energies. For realization 21, 
which has two well separated magnetization peaks, there are also two well separated peaks 
in the bond energy distribution. Similarly, for realizations 25 and 31 there are two peaks 
in the bond energy distribution. In a first-order interpretation of these realizations the + 
and — phases have the lower energy and the phase the higher. On the other hand, for 
realization 14, which has five magnetization peaks, there is a single broad maximum in the 
bond energy histogram with a significant, low energy shoulder. Altogether, 5 out of 22 
realizations show two well defined peaks in the bond energy histogram while the remaining 
realizations show a single peak though sometimes with a significant low energy shoulder. 

The joint magnetization/bond-energy distribution for realization 25 is shown in Fig. |15. 
For comparison, the joint distribution for the pure Ising model at the infinite system critical 
temperature is shown in Fig. |16|. The three lobes in the distribution correspond to the three 
peaks in the magnetization histogram. Fig. ^ It is tempting to interpret realizations 25 and 
31 in the language of first-order transitions and declare the low energy side lobes as ordered 
phases and the central lobe as a disordered phase. 

The importance of fine-tuning the magnetization and temperature for each realization 



of disorder is illustrated in Figs. [T^, 0, |19| and The magnetization histogram for real- 



ization 1 at the fine-tuned parameter values, Pc = 0.268385 and He = 0.00127049 is shown 
in Fig. 17. Note, that this value of He is not the same as the value that would be ob- 
tained by summing the local fields. For realization one, that value is almost 50% larger, 
—A J2i hi/N = 0.00178774. Figure |18| shows realization 1 with same external field but with 
the temperature raised by 5% leaving only the disordered phases. Figure O shows realiza- 
tion 1 with the same external field but the temperature decreased by 5%. Although, the 
external field is too negative here, it is nonetheless clear that the two ordered phases are 
now dominant over the disordered phase, which has nearly vanished. Finally, if the external 
field is set to and the feedback temperature is used, only the negative peak survives, as 



shown in Fig. EO 



We considered several realizations at stronger disorder. Figure ^ shows the magne- 
tization histogram for realization 1 at disorder A = 0.433 at the value of external field 
and temperature determined from the fine-tuning procedure. Comparing Fig. ^ to Fig. |17 
we see that the three-peaked structure is conserved and the peaks becomes sharper as the 
strength of the disorder is increased. In this experiment, we used 16 replicas with maxi- 
mum disorder strength, A = 0.5 however, the fine-tuning procedure failed for the two most 
disordered replicas and only yielded the minus phase. In general, the feedback procedure 
was less stable for stronger disorder. Figure ^ shows realization 14 at disorder A = 0.5 at 
the values of the temperature and external field determined from the fine-tuning procedure. 
A comparison to Fig. |^ shows that increasing the disorder has sharpened the five peaked 
structure of the magnetization histogram. It would be interesting to carry out a systematic 
study of stronger disorder to see if there is a trend toward more "first-order" like realization 
however this will require developing a better fine-tuning mechanism. 



7 



Overall, out of 21 realizations (seeds 10-30) at disorder strength A = .35 we found that 
the average values of the critical parameters is {{3^ = 0.2663 and (He) = 0.0006 (compared 
with initial values before the feedback procedure of /3 = 0.2670 and H = 0) with a standard 
deviations of = 0.0037 and an = 0.0031. The standard deviation of He is consistent with 
the anticipated value an ~ AL"'^/^ ^ 0.0030. 



B. Dynamics of the algorithm 



We have studied the dynamics of the replica exchange algorithm and compared it to the 
Metropolis algorithm. There are two ways of comparing the algorithms. First, time can 
be measured by Monte Carlo sweeps. This approach favors the replica exchange algorithm 
since a single Monte Carlo sweep involves flipping spins in many (here 16) replicas and since 
growing clusters is computationally more intensive than flipping single spins. The second 
approach is to compare actual running time for the two algorithms on the same computer. 
This approach ignores the possibility that the one algorithm is better optimized and it is 
not an appropriate way to measure fundamental quantities such as the dynamic exponent 
however it does give a practical comparison and is useful for deciding which algorithm to 
choose for a given system size. For 24^ systems with 16 replicas with maximum disorder 
A = .35 on a 450MHz Pentium III machine we find that the replica exchange algorithm runs 
at about 0.56 sec/sweep and that our implementation of the Metropolis algorithm runs at 
about 0.0088 sec/sweep. As is the case for the equilibrium properties of the system, we find 
great differences in the dynamics depending on the realization of disorder. These differences 
are much more pronounced for the Metropolis algorithm. We measured the integrated 
autocorrelation time for the magnetization for realization 20 for both algorithms with the 
result r^^*™P°"" = 6000 and t'^^^'^^ = 1100 when measured in Monte Carlo sweeps. The 
advantage of the replica exchange algorithm is lost however when these times are converted 
to running times; T-^ietropohs _ ^rephca _ QiQ^ fjj^g joint energy/magnetization 

histogram for realization 20 is smooth and without gaps suggesting that there are no global 
energy barriers separating regions of phase space. At the other extreme is realization 31 for 
which the joint energy/magnetization distribution is split into three pieces separated by gaps 
where the probability is very small. Using the parameters for realization 31 determined by 
fine-tuning for A = .3267, we did a long Metropolis run of 100 million Monte Carlo sweeps 
starting from random initial spins. The simulation stayed mainly in the "0" state with two 
excursions to the "+" state. The "— " state was never visited. On the other hand, two 
Metropolis runs of 20 million sweeps starting from all spins up and all spins down went to 
the "0" and "— " states, respectively, and stayed there for the entire simulation. Figure ^ 
shows the magnetization time series from the replica exchange algorithm for realization 31 
for 800,000 sweeps. Although it is difficult to estimate the autocorrelation time from this 
series, it is clear that the simulation samples each of the states many times. The conclusion 
is that for most realizations of disorder, at size 24^ there is no great practical advantage 
to using the replica exchange algorithm but that for cases where well separated regions 
are present in the joint distribution, only the replica exchange algorithm is able to reach 
equilibrium within reasonable simulation times. 
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IV. DISCUSSION 



One of the chief motivations for this study was to determine the nature of the RFIM 
phase transition. Unfortunately, our resuhs fail to settle this question. What would we 
have expected from the different scenarios for the phase transition? If the transition is 
first-order we should find three phases in coexistence at the transition point. The high 
temperature phase, would have little magnetization and a large bond energy while the two 
low temperature phases would have large absolute values of magnetization, one positive and 
the other negative, and small bond energies. This situation is exactly what is seen in roughly 
10% of realizations. For systems much smaller than 24^, three phase coexistence is not seen 
simply because the magnetization distribution is too broad to display three distinct peaks. 
One hypothesis is that, as system size becomes large, the fraction of systems displaying 
three phases increases and that, in the thermodynamic limit, the transition at A = .35 is 
first-order. However, the existence of realization 14 with its five peaks is difficult to reconcile 
with this viewpoint. 

If the transition is continuous and follows the droplet model scenario |l6|,|8|j25[| we would 
expect two peaks in the magnetization histogram corresponding to two states of a single 
critical phase. The width of each peak should scale as L~"^l^ and the separation between 
the peaks should scale as L~^l^ . The majority of our disorder realizations display two peaks 
and an alternative to the first-order transition hypothesis is that, as system size increases, 
the fraction of realizations with two peaks in the magnetization histogram approaches 100%. 
Although the original droplet model envisioned two states with roughly the same energy and 
differing by flipping a critical cluster of spins, it is possible to imagine more complicated 
droplet pictures where there are sometimes more than two states in the critical phase and 
therefore more than two peaks in the magnetization histogram. However, our five realizations 
that display more than one peak in the bond energy histogram are difficult to reconcile 
with the droplet model since one expects the states, however many there are, to be nearly 
degenerate in bond energy. 

Ultimately, the essential difference between the first-order scenario and the droplet model 
continuous-transition scenario is whether there are distinct phases at the transition point 
or whether the peaks in the magnetization histogram are different states that are part of a 
single critical phase. The comes down to the question of whether the peaks in magnetization 
histogram move toward the origin as L — >■ cxd or remain at finite values. Unfortunately, the 
consensus is that the magnetization exponent, if it is not actually zero, is very tiny so 
that the decrease in magnetization with system size would be far too weak to observe in 
simulations. 



V. CONCLUSIONS 

We have studied the phase transition of the random field Ising model by Monte Carlo 
simulation of systems of size 24^ using a new cluster algorithm. We find major qualitative 
differences between the behavior of different realizations of disorder. Some realizations 
display a broad two-peaked magnetization histogram consistent with a continuous transition, 
while a small fraction display a three-peaked structure consistent with a first-order transition. 
Our main conclusion is that more work needs to be done to determine the nature of the 
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transition. It would be very useful to study larger system sizes and stronger disorder to 
determine whether there are trends in the qualitative features of the phase transition. Does 
the fraction of "first-order" systems increase as disorder or system size increases? It will 
require improvements of the algorithm including the fine-tuning technique or substantially 
more computer power to resolve these questions. 
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FIGURES 
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FIG. 1. Proposed phase diagram for the 3D RFIM in the H-T plane at a fixed strength of 
disorder. The bold lines are first-order lines. The black dot is the thermal first-order transition at 
{Tc,H = 0) where the -|-,— and phases coexist. The open circles are critical endpoints of the two 
order-disorder first-order lines. 




13 



FIG. 2. Approximate phase diagram for the RFIM and rephcas for the rephca exchange method. 
The phase boundary in the 1//5-A//? plane is taken to be an chiptical curve starting at the 
pure Ising critical point {1/(3 = 4.512, A//3 = 0) and ending at the zero temperature transition 
(1//3 = 0, A//3 = 2.35). The initial conditions for the 16 replicas lie on this curve and are evenly 
space in A. During the feedback process, each replica is shifted by a small amount in (3 and H. 
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FIG. 4. Magnetization histogram for realization 21 at A = 0.35. 
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FIG. 5. Magnetization histogram for realization 25 and A = 0.35. 
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FIG. 6. Magnetization histogram for reahzation 31 at A = 0.3267. 
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FIG. 8. Magnetization histogram for the pure Ising model at the infinite system size critical 
temperature Pc = 0.22165. 
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FIG. 11. Bond energy histogram for realization 25 at A = 0.35. 
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FIG. 14. Bond energy histogram for the pure Ising model at the infinite system size critical 
temperature (5c = 0.22165. 
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FIG. 15. The joint magnetization/bond-energy distribution for realization 25 at A = 0.35. 
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FIG. 16. The joint magnctization/bond-cnergy distribution for for the pure Ising model at the 
infinite system size critical temperature (3c = 0.22165.. 
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FIG. 17. Magnetization histogram for realization 1 at A = 0.35 and the phase transition point 
determined from fine-tuning, /3c = 0.268385 and = 0.00127049. 





60 




50 




40 






1 


30 








20 




10 









1 



1 



-0.5 0.5 

magnetization per spin 

FIG. 18. Magnetization histogram for reaHzation 1 at A = 0.35 and /3 = 0.267041 (5% warmer 
than Pc determined for reahzation 1) and He = 0.00127049. 



21 



100 



42 



80 
60 
40 
20 





-1 -0.5 0.5 

magnetization per spin 

FIG. 19. Magnetization histogram for realization 1 at A = 0.35 and /? 
than Pc determined for reahzation 1) and He = 0.00127049. 
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FIG. 20. Magnetization histogram for realization 1 at A = 0.35, (3c = 0.267041 and H 
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FIG. 21. Magnetization histogram for realization 1 at A = 0.433 at temperature and external 
field determined from fine-tuning, /3c = 0.289176 and He = 0.00158005. 
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FIG. 22. Magnetization histogram for realization 14 at A = 0.5 at temperature and external 
field determined from fine-tuning. 
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FIG. 23. Magnetization time series for realization 31 for A = .3267. 
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